MTA Ridership & Turnstile‑Hopping Data Visualization¶

# Load data from CSV file into a pandas DataFrame
df = pd.read_csv('clean_subway.csv')

df.sample(10)

print(df['station_complex_id'].value_counts())
df[df['station_complex_id'] == 76]

len(df)

df.drop_duplicates(inplace=True)
len(df)

# Create an hourly timestamp column as integer
df['hour'] = df['transit_timestamp'].dt.floor('H')

# Group by station and hour, calculating mean hourly ridership
hourly_data = df.groupby(['station_complex', 'hour', 'latitude', 'longitude'])['ridership'].mean().reset_index()

# Calculate percentiles for normalization
low, high = np.percentile(hourly_data['ridership'], [5, 95])

# Create figure and axis
fig, ax = plt.subplots(figsize=(12, 8))

# Set up the plot limits based on your data
ax.set_xlim(df['longitude'].min(), df['longitude'].max())
ax.set_ylim(df['latitude'].min(), df['latitude'].max())

# Add a background map using contextily
ctx.add_basemap(ax, crs="EPSG:4326", source=ctx.providers.CartoDB.Positron)

# Initialize scatter plot with RdYlGn_r colormap
scatter = ax.scatter([], [], c=[], cmap='RdYlGn_r', alpha=0.6, s=50, vmin=low, vmax=high)

# Set up the plot
ax.set_title('NYC Subway Hourly Ridership in 2024', fontsize=16)
ax.set_xlabel('Longitude', fontsize=12)
ax.set_ylabel('Latitude', fontsize=12)

# Add colorbar
cbar = plt.colorbar(scatter)
cbar.set_label('Average Hourly Ridership', fontsize=12)

# Add timestamp text
timestamp_text = ax.text(0.02, 0.98, '', transform=ax.transAxes, fontsize=12, verticalalignment='top', 
                         bbox=dict(facecolor='white', alpha=0.7, edgecolor='none'))

# Get unique hours for animation
unique_hours = sorted(hourly_data['hour'].unique())

# Animation update function
def update(frame):
    current_hour = unique_hours[frame]
    current_data = hourly_data[hourly_data['hour'] == current_hour]
    
    scatter.set_offsets(current_data[['longitude', 'latitude']])
    scatter.set_array(current_data['ridership'])
    
    timestamp_text.set_text(current_hour.strftime('%Y-%m-%d %H:%M'))
    
    return scatter, timestamp_text

# Create animation
ani = animation.FuncAnimation(fig, update, frames=len(unique_hours), 
                              interval=200, blit=True)

plt.tight_layout()
plt.show()

# Uncomment the following line to save the animation as gif (requires pillow)
ani.save('nyc_subway_hourly_ridership_in_2024_with_map.gif', writer='pillow', fps=20)
print("Animation saved successfully.")

# Create an hourly timestamp column as datetime
df['hour'] = df['transit_timestamp'].dt.hour

# Group by station and hour, calculating mean hourly ridership
hourly_data = df.groupby(['station_complex', 'hour', 'latitude', 'longitude'])['ridership'].mean().reset_index()

# Calculate percentiles for normalization
low, high = np.percentile(hourly_data['ridership'], [5, 95])

# Create figure and axis
fig, ax = plt.subplots(figsize=(12, 8))

# Set up the plot limits based on your data
ax.set_xlim(df['longitude'].min(), df['longitude'].max())
ax.set_ylim(df['latitude'].min(), df['latitude'].max())
ax.set_title('NYC Subway Average Hourly Ridership', fontsize=16)
ax.set_xlabel('Longitude', fontsize=12)
ax.set_ylabel('Latitude', fontsize=12)

# Add a background map using contextily with OpenStreetMap
ctx.add_basemap(ax, crs="EPSG:4326", source=ctx.providers.OpenStreetMap.Mapnik)

# Initialize scatter plot
scatter = ax.scatter([], [], c=[], cmap=cmap, alpha=0.6, s=50, vmin=low, vmax=high)

# Add colorbar
cbar = plt.colorbar(scatter)
cbar.set_label('Average Hourly Ridership', fontsize=12)

# Add hour text
hour_text = ax.text(0.02, 0.98, '', transform=ax.transAxes, fontsize=12, verticalalignment='top')

# Animation update function
def update(frame):
    current_hour = frame
    current_data = hourly_data[hourly_data['hour'] == current_hour]
    
    scatter.set_offsets(current_data[['longitude', 'latitude']])
    scatter.set_array(current_data['ridership'])
    
    hour_text.set_text(f'Hour: {current_hour:02d}:00')
    
    return scatter, hour_text

# Create animation
ani = animation.FuncAnimation(fig, update, frames=24, interval=500, blit=True)

plt.tight_layout()
plt.show()

# Save animation
ani.save('nyc_subway_average_hourly_ridership_with_map.gif', writer='pillow', fps=2)
print("Animation saved successfully.")

display(Image(filename='nyc_subway_average_hourly_ridership_with_map.gif'))

# Create an hourly timestamp column as datetime
df['hour'] = df['transit_timestamp'].dt.hour

# Group by station and hour, calculating mean hourly ridership
hourly_data = df.groupby(['station_complex', 'hour', 'latitude', 'longitude'])['ridership'].mean().reset_index()

# Calculate percentiles for normalization
low, high = np.percentile(hourly_data['ridership'], [5, 95])

# Determine fixed limits for the bar chart based on overall data
overall_min_ridership = hourly_data['ridership'].min()
overall_max_ridership = hourly_data['ridership'].max()

# Create figure and axes for scatter plot and bar chart
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(18, 8))

# Set up the plot limits based on your data for scatter plot
ax1.set_xlim(df['longitude'].min(), df['longitude'].max())
ax1.set_ylim(df['latitude'].min(), df['latitude'].max())
ax1.set_title('NYC Subway Average Hourly Ridership', fontsize=16)
ax1.set_xlabel('Longitude', fontsize=12)
ax1.set_ylabel('Latitude', fontsize=12)

# Add a background map using contextily with OpenStreetMap
ctx.add_basemap(ax1, crs="EPSG:4326", source=ctx.providers.OpenStreetMap.Mapnik)

# Initialize scatter plot
scatter = ax1.scatter([], [], c=[], cmap=cmap, alpha=0.6, s=50, vmin=low, vmax=high)

# Add colorbar to scatter plot
cbar = plt.colorbar(scatter, ax=ax1)
cbar.set_label('Average Hourly Ridership', fontsize=12)

# Add hour text to scatter plot
hour_text = ax1.text(0.02, 0.98, '', transform=ax1.transAxes, fontsize=12, verticalalignment='top')

# Animation update function
def update(frame):
    current_hour = frame
    current_data = hourly_data[hourly_data['hour'] == current_hour]
    
    # Update scatter plot data
    scatter.set_offsets(current_data[['longitude', 'latitude']])
    scatter.set_array(current_data['ridership'])
    
    # Update hour text
    hour_text.set_text(f'Hour: {current_hour:02d}:00')
    
    # Get top 10 stations by ridership for the current hour
    top_stations = current_data.nlargest(10, 'ridership')
    
    # Clear previous bar chart and plot new data
    ax2.clear()
    ax2.barh(top_stations['station_complex'], top_stations['ridership'], color='skyblue')
    ax2.set_title(f'Top 10 Stations - Hour: {current_hour:02d}:00', fontsize=16)
    ax2.set_xlabel('Average Hourly Ridership', fontsize=12)
    
    # Set consistent x-axis limits for the bar chart based on overall data
    ax2.set_xlim(overall_min_ridership * 0.9, overall_max_ridership * 1.1)  # Adjust margins as needed
    
    ax2.invert_yaxis()  # Invert y-axis to have the highest ridership on top
    
    return scatter, hour_text

# Create animation
ani = animation.FuncAnimation(fig, update, frames=24, interval=500, blit=True)

plt.tight_layout()
plt.show()

# Save animation
ani.save('nyc_subway_average_hourly_ridership_with_map_and_graph.gif', writer='pillow', fps=2)
print("Animation saved successfully.")

display(Image(filename='nyc_subway_average_hourly_ridership_with_map_and_graph.gif'))

arrests_df = pd.read_csv('NYPD_Arrests_small.csv')
subway_df = pd.read_csv('clean_subway.csv')

# Create a new dataframe with only the columns that have THEFT OF SERVICES in the PD_DESC column
turnstile_jumps = arrests_df[arrests_df['PD_DESC'].str.contains('THEFT OF SERVICES', case=False, na=False)]


len(turnstile_jumps)
turnstile_jumps.sample(5)

def find_nearest_station(lat, lon, stations):
    arrest_coords = np.array([[lat, lon]])
    station_coords = stations[['latitude', 'longitude']].values
    distances = cdist(arrest_coords, station_coords)
    nearest_index = distances.argmin()
    return stations.iloc[nearest_index]['station_complex']

# Apply function to get nearest station for each arrest
turnstile_jumps['nearest_station'] = turnstile_jumps.apply(
    lambda row: find_nearest_station(row['Latitude'], row['Longitude'], subway_df),
    axis=1
    )

# Count turnstile jumpings per station
jumps_per_station = turnstile_jumps['nearest_station'].value_counts().reset_index()
jumps_per_station.columns = ['station_complex', 'turnstile_jumpings']

# Merge with subway data to get coordinates
result_df = pd.merge(jumps_per_station, subway_df[['station_complex', 'latitude', 'longitude']], on='station_complex', how='left')
result_df = result_df.drop_duplicates(subset=['station_complex'], ignore_index=True)

# Select relevant columns for the enriched arrests data
enriched_arrests = turnstile_jumps[['ARREST_KEY', 'ARREST_DATE', 'PD_DESC', 'ARREST_BORO', 'Latitude', 'Longitude', 'nearest_station']]

# Export enriched arrests data
enriched_arrests.to_csv('enriched_turnstile_jumping_arrests.csv', index=False)
print("Exported 'enriched_turnstile_jumping_arrests.csv'")

# Dataset only goes from October to December 2023
enriched_arrests
enriched_arrests['ARREST_DATE'].agg(['min', 'max'])

result_df

display(nyc_map)